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Abstract 

Analytical models have been developed for fracture propagation over the last several 
decades and are now considered with renewed interest; the range of their applicability varies 
for different materials and different loading conditions. Systematic experimental measure¬ 
ments and numerical simulations are presented against an asymptotic prediction for the 
crack trajectory developing in solids with small voids (of circular and elliptical shape). The 
experiments also cover the dynamic regime, where new features involving crack kinking and 
roughness of the fracture surface occur. 
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1 Introduction 

Fracture is the most limiting factor for the development of highly brittle materials such as ad¬ 
vanced ceramics. For these materials, theoretical predictions of the geometrical shape of cracks 
are rare and when such predictions are possible they are always valuable both for design pur¬ 
poses and for the advances in the knowledge of these fragile media. The deflection of a crack 
from a straight line is naturally observed in experiments and linked to an increase of toughness 
laiisiiiH]. Crack trajectories may be deflected from rectilinearity due to an interaction with 
imperfections such as small voids or inclusions, which may result from exploitation damage or 
manufacturing defects. Crack branching and path deflection in two-dimensional elastic solids 
were examined by Sumi et al. (1983) [21], where the asymptotic method was also illustrated 
for the case of a slightly slanted Griffith crack under bi-axial load. Computational models have 
been developed by Xu et al. (1994; 1998) [28][29] to study the crack growth in heterogeneous 
solids, while an energy minimizer technique for crack path determination has been developed 
by Francfort and Marigo (1998) [7] and Dal Maso et al. (2005) [6]. The local higher-order 
asymptotic analysis for the elastic fields near the crack tip was developed by Sumi (1986) [23] 
and the idea of using the higher-order terms to predict a fracture path was implemented in 
the computational FEM algorithm [22]: to our knowledge this was the first systematic attempt 
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towards the formulation of the hybrid method, based on superposition of analytical and numer¬ 
ical solutions, employing the coefficients in the hrst-order perturbation formulae. The paper 
by Sumi (1986) [23] presented simulations for both kinking and a smooth crack advance. The 
hnite element simulation for a system of cracks, based on the hrst-order perturbation method 
and an incremental quasi-static crack growth assumption, was developed by Sumi and Wang 
(1998) [25]. Formulations for a curved crack based on a perturbation procedure were discussed 
by Hori and Vaikuntan (1997) [H], although no weight functions were involved in their analysis. 

Two-dimensional models (both plane stress and plane strain) of crack propagation have been 
developed in the analytical form, in which the crack path was predicted in terms of elementary 
functions or via a straightforward integration involving weight functions. These approaches 
include an asymptotic model describing the interaction of a semi-inhnite crack with small voids 
[I5],[IS!,[I3],[I1!, where the voids are characterized by their dipole tensors m, also known as 
the virtual mass tensor. In the case of vector elasticity, the dipole tensors are of rank 4, but 
allow for a matrix representation in a specially chosen basis; the corresponding matrices are 
referred to as Polya-Szego matrices mi. 

Practical needs require efficient tools, and the asymptotic model, appealing for its simplicity, 
can be used for many fracture conhgurations. In fact, Bigoni et al. (1998) [2] and Valentini et 
al. (1999, 2002) [27] [26] have successfully compared theoretical predictions with experimental 
measurements, but in an indirect way, i.e. the cracks were induced by an indenter, which 
provided a complex stress state, different from the remote Mode-I condition assumed in the 
model. In the present article, the reliability of the asymptotic model is assessed for quasi-static 
regimes, together with the observation of important effects related to dynamic crack kinking. 
In particular, in addition to the asymptotic model, a computational technique, namely, the 
extended hnite element method (X-FEM) as implemented in Abaqus FEA (Shi et al., 2010 
m), is employed to predict crack propagation. This method (already used for the analysis of 
structural components failure m, 0 and of quasi-static crack growth m, mi) is employed, 
together with the asymptotic model, to predict crack paths in brittle (PMMA) notched plates 
subject to Mode-I loading and containing elliptical voids. An excellent agreement between 
theoretical results and experiments is observed, with the asymptotic solution correctly predicting 
fracture dehection toward the small void, as shown for instance in Fig. [TJ where a fracture path 
deviates from the straight line due to the interaction with thin elliptical voids. In the hgure 
the cyan/dashed curve represents the shape predicted by the asymptotic model, whereas the 
solid/red curve shows the actual shape obtained in the experiment (the FE simulation coincides 
with the asymptotics and hence is not reported in the hgure). 

From this result it can be deduced that the asymptotic approach fully explains crack de¬ 
hection induced by a void (when sufficiently far from the fracture trajectory to satisfy the 
asymptotic assumptions), so that it can be employed in design situations involving the presence 
of defects in ceramics. 

2 Modelling 

2.1 Asymptotic model 

In this section we give a brief outline of the model used for the analysis of crack trajectories 
inhuenced by elliptical voids; for a more detailed description see Movchan et al. (1991)[IS], 
Movchan (1992)[T2], Movchan and Movchan (1995) [U], Movchan et al. (1997)[IH]. Valentini et 
al. (1999)123. 

An infinite, brittle isotropic linearly elastic body is assumed, containing a semi-infinite crack 
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Figure 1: Experimental crack trajectory (a post-mortem photo of a PMMA sample loaded in Mode-I) compared 
to the theoretical prediction (reported with a cyan/dashed line) for two cases of interaction with one (A) and two 
(B) elliptical voids. The crack trajectory is highlighted by means of a solid/red colored strip inserted between 
the two crack surfaces. The parameters indicated in the figure are explained in Fig. [2| 


interacting with some defects (voids or inclusions), and propagating under pure Mode-I loading, 
corresponding to a stress intensity factor Kj greater than the critical one. The elastic properties 
are specihed through the Lame constants A and fi for the body and Aq and /io for the inclusion. 
We assume that the defects present in the inhnite elastic medium are Tar’ from the straight 
trajectory that would be followed by the crack in the absence of disturbances (Fig.|^. The 




Figure 2: The analyzed crack geometry interacting with an elliptical void (A) and an elliptical elastic inclusion 
(B). 


shape of the defect is selected to be an ellipse centered at the point (x^, x^) with the major and 
minor semi-axes denoted by a and 6, respectively. The major axis is inclined at an angle 9 with 
respect to the xi-axis. Due to the presence of the defect the crack trajectory may deviate from 
a straight line and this perturbation can be described by the closed-form formula for the crack 
trajectory i?(/), expressed as a function of the crack tip coordinate I (Valentini et ah, 1999; 
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The elastic properties of the inclusion and matrix influence the crack trajectory via the 
constants: 
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where k. — 3 — Au (for plane strain) or = (3 — z/)/(l + z/) (for plane stress), and z/ G ( — 1,1/2) 
is the Poisson’s ratio. 

In the case of plane stress^ the formula for the crack trajectory due to elliptical void may be 
obtained by taking /io = Aq = 0 in 0-@> which gives 
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and m = - 

a + b 


It is important to remark that expression Q depends on the Poisson’s ratio z/, the angle of 
inclination 9 of the major axis, and the parameters R and m. 

We note that imperfectly bonded circular inclusions have been considered in Bigoni et al. 
(1998)[2]. 


2.2 Adaptive numerical simulation 

An algorithm has been implemented in ABAQUS to simulate a quasi-static crack growth in 
a solid with small voids. It enables a crack advance, according to a brittle fracture criterion, 
without any remeshing. 


4 























Published in International Journal of Solids and Structures (2015) 63 219-225 

doi:10.1016/j.ijsolstr.2015.03.001 


We have computed the trajectory for crack propagation by means of the 2D extended fi¬ 
nite element method (X-FEM), implemented in the Abaqus/Standard hnite element package. 
The extended Finite Element Method is a partition of unity based method, particularly suit¬ 
able for modelling crack propagation phenomena, without a-priori knowledge of the crack path 
and without remeshing. The simulations were performed using a parametric python script for 
ABAQUS, that is run by means of MatLab. CPS4 elements were used for 2D plane stress 
analyses. 



Figure 3: X-FEM Abaqus simulations under the plane stress assumption for one (A) and two (B) circular voids 
are reported (solid/black line) and compared to results obtained with the asymptotic formulation (cyan/dashed 
line). 


The crack trajectory evolution was modelled via Linear Elastic Fracture Mechanics approach 
(LEFM) specihc for brittle fracture. In particular, the strain energy release rate at the crack tip 
was calculated using the virtual crack closure technique (VCCT), so that the crack propagates 
when the strain energy release rate reaches a critical value. We have chosen the so-called 
‘BK law’ to simulate mix-mode behavior (tolerance=0.1, viscosity=0.0001) and the Maximum 
Tangential Stress (MTS) as the criterion for singling out the crack propagation direction. The 
former criterion is a combination of the energy release rates in Mode I, Mode II, and Mode 
III in a power law relationship (Benzeggagh and Kenane, 1996) [I]. The latter criterion uses 
the maximum tangential stress to determine the normal direction to the crack plane (Erdogan 
and Sih, 1963) |1]. Note that the X-FEM approach allows us to avoid the use of special hnite 
elements with embedded singularities, which requires that the crack propagates across an entire 
element at a time. The strain energy release rate at the crack tip has been computed by the 
modihed Virtual Crack Closure Technique (VCCT) in combination with principles of linear 
elastic fracture mechanics. The VCCT is based on the assumption that the strain energy 
released for crack extension of a certain amount is the same as the energy required to close the 
crack of the same amount. 

The computation set-up, in term of sample size, the way to apply the load and the material 
properties were selected to be identical to the experimental set-up; in particular, the load is 
applied by imposing a displacement at the grips of 0.8 /im/s, in order to maintain a low loading 
rate. 

Two examples of the X-FEM crack path predictions are shown in Fig. [^and compared to 
the prediction of the asymptotic model, eq. Q, for one and two circular voids. An excellent 
agreement can be observed. 
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An asymptotic formula for the crack trajectory interacting with a small void, characterized 
by a 3 X 3 dipole matrix V is given in m, see also US]. The following formula is valid for a 
crack interacting with an arbitrary small void (not necessarily elliptical) 


4 /i 
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In the specific case of an elliptical void, the dipole matrix V becomes (Movchan et al. 

(1997)[IS|) 


p(«)= 


/(K-l)[S-i?2(K-l)]-S 

_ 




^ \ 
K — 1 


(K-ir 

{k-1)[-E-R^{k-1)]-E a 




(n — 1)^ 
A 


n — 1 


A 

n — 1 


K, — 1 

-2i?2 


(7) 


where 


q = 


A + /i 


S = 4mi?^ cos20, 


( 8 ) 
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The crack deflection at inflnity, Tar away’ from the small void, can be obtained from eq. (§ 
in the limit / ^ oc, which for a generic void becomes 


H{co) = ^[£*(0)iP£(0)], 


and in the particular case of an ellipse 


K + 1 


R 


H{oo) = — 0(1 + rn ). 


(9) 


( 10 ) 


We would like to emphasize that, although the deflection at inflnity depends on the dipole 
matrix P, it is independent of the void’s orientation (independent of 9 for the elliptical void). 
This important feature is clearly visible from both the Abaqus simulations and the experiments 
(see Figs. mil and[^. 


3 Experiments 

The experimental veriflcation of the theoretical asymptotic model outlined above was performed 
(at the ‘Instabilities Lab’, http://www.ing.unitn.it/dims/ssmg/) by quasi-statically loading 
(controlled displacements) notched plates of a brittle material under Mode-I conditions un¬ 
til failure. The post-mortem crack trajectories were analyzed and compared to the predictions 
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of the asymptotic method and to those of the numerical model. We highlight the limitations of 
the asymptotic approach to the cases of quasi-static crack advance, showing that in the dynamic 
regime an instability may occur, leading to a crack kinking. The conhgurations reported in the 
text include several geometries of small voids at different orientations. The implementation of 
the experiment was done to ensure the required quasi-static growth of a crack. 

3.1 Quasi-static crack propagation 

PMMA (PERPEX white and black from Bayer, with elastic modulus 3350 MPa, z/=0.38) plates 
with a V-shaped notch (to be loaded in-plane) were employed for the study of crack propaga¬ 
tion, in agreement with the standard test method ASTM E647-00. Samples of two different 
dimensions were used, namely, 125 mm x 95.0 mm x 3.0 mm and 130 mm x 105 mm x 3mm, for 
testing one or more circular or elliptical voids, respectively. The notches have been engraved 
with a cutter (realized from a blade of a Stanley 99E utility knife to advance of a hxed amount 
of 0.4 mm with a micrometrical screw, see Eig. @0) to trigger the propagation of the fracture. 
Specimens were obtained by cutting a PMMA plate with a EGX-600 Engraving Machine (by 
Roland). The position and the size of the voids are reported in the labels of the hgures, with 
reference to the nomenclature introduced in Eig. 0A. 

Tensile force on the structure was provided by imposing displacement with a load frame MIDI 
10 (Messphysik) and the load measured with a load cell Gefran TH-KN2D. Data was acquired 
with a NI CompactDAQ system, interfaced with Labview 8.5.1 (National Instruments). Photos 
were taken with a Sony NEX 5N digital camera (equipped with 3.5-5.6/18-55 lens, optical steady 
shot from Sony Corporation), with a Nikon D200 digital camera (equipped with a AE-S micro 
Nikkor lens 105 mm 1:2.8G ED), or with a stereoscopic microscope Nikon SMZ-800 (equipped 
with objectives P-Plan Apo 0.5X and P-ED Plan 1.5X). A very low rate of loading, namely, 0.8 
fim/s was imposed during the test to ensure quasi-static crack propagation. 

3.2 Experimental measurements versus analytical prediction 

Eor the samples where small (circular and elliptical) voids he at different relative positions 
and orientations, but always sufficiently distant, from the propagating crack, we demonstrate 
an impressive agreement between the theoretical prediction of the asymptotic model and the 
experimental measurements. It is also noted that the crack deflection at inhnity is approximately 
the same for different orientations of the voids of non-circular geometry, which is in agreement 
with the prediction of the asymptotic model. 

3.2.1 Circular versus elliptical void attracting the crack 

Experiments for both circular and elliptical voids, were performed, Eigs. aim and compared 
to theoretical results obtained with the asymptotic approach, eq. 

The numerical simulations were found to match the the predictions of the asymtotic so well 
that they were superimposed on the predictions of the asymptotic model and therefore have 
not been reported. 

Two experiments for one and two voids are reported in each of Eigs. [^an di elliptical voids 
were considered in the former hgure and circular in the latter. Beside the very satisfactory 
theoretical/experimental agreement, it is important to remark that, in the case of two voids 
(elliptical or circular), the deflection at inhnity is close to zero, as a consequence of the sum of 
two trajectory dehections with opposite sign. 
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A) 


Figure 4: Experimental crack trajectory (a post-mortem photo of a PMMA sample loaded in Mode-I) compared 
to the theoretical prediction (reported with a cyan/dashed line) for two cases of interaction with one (A) and 
two (B) circular voids. The crack trajectory is highlighted by means of a solid/red colored strip inserted between 
the two crack surfaces. (C) Detail of the engrave cut at the notch tip to trigger rectilinear fracture propagation. 




The influence of the orientation of an elliptical void has been investigated in Fig. where 
three experiments are reported, all confirming the theoretical prediction that the inclination of 
the ellipses does not influence crack deflection at infinity. 

In conclusion, the two most important phenomena that are very evident from the experi¬ 
ments are that the voids ‘attract’ the crack trajectory and that the deflection at infinity only 
depends on the dipole matrix as predicted from the asymptotic model, eq. (10). 


3.3 Applicability of the asymptotic model, dynamics and crack kinking 

In this section we highlight the limitations of the above-described asymptotic model in the case 
when the crack trajectory is too close to the void and when propagation of the crack reaches 
the dynamic regime. 

The case of a large circular void close to the rectilinear crack trajectory (that would occur in 
the absence of the void) was investigated and results are reported in Fig. for the experimental 
fracture trajectory (Fig. A) and its asymptotic and numerical predictions (Fig. [^B). The 
asymptotic model is based on the assumptions that the ratio e between the inclusion diameter 
and the distance from the inclusion centre to the crack trajectory is small, an assumption now 
clearly violated. It is interesting to notice that the crack trajectory can be captured with the 
numerical model, while the asymptotic evaluation still provides a good approximation to the 
initial crack deviation. 

The other departure from the hypotheses of both the asymptotic and the numerical approach 
occurs when a sufficiently high (and constant) rate of loading is applied to the sample. In this 
case, an acceleration of the crack occurs, leading subsequently to a dynamic kinking of the 
crack trajectory. The latter differs from the case of the quasi-static crack propagation, where 
the crack path remains smooth; the loading rate has to be reduced significantly (as done in 
the previously presented experiments) in order to achieve conditions appropriate for a quasi¬ 
static crack growth. From the experiments we have seen that there is a wide range of applied 
displacement rate within which the asymptotic model gives a good prediction of the quasi¬ 
static crack trajectory. In particular, if the loading rate determined by the speed of the loading 
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Figure 5: Experimental crack trajectory (a post-mortem photo of a PMMA sample loaded in Mode-I) compared 
to the theoretical prediction (reported with a white cyan/dashed line) for two cases of interaction with one 
elliptical void: major axis is inclined at an angle ^ = 0° (A), 0 = 90° (B) and 0 = 45° (C) with respect to the 
xi-axis. The crack trajectory is highlighted by means of a solid/red colored strip inserted between the two crack 
surfaces. 


movable crosshead (shown in Fig. is constant and exceeds 10 mm/s, then the crack growth 
accelerates and becomes unstable, resulting in a roughness of the fracture surface and kinking of 
the crack trajectory. On the other hand, if the loading rate falls below 10 mm/s, the crack path 
remains smooth and follows the theoretical prediction given by the asymptotic and numerical 
models, as previously shown. Fig. shows kinking instability in the case of dynamic crack 
growth as related to the roughness of the crack path shown in Fig. (B), to be compared to 
the smoothness of the paths found in ah cases of quasi-static propagation (see for example Fig. 
[|A). 

4 Conclusions 

An experimental and computational (using adaptive hnite elements with an X-FEM algorithm 
and moving boundary associated with an advancing crack) verihcation has been presented 
for an asymptotic model, developed to predict the path of a crack propagating in a brittle 
material. In the quasi-static regime, the predictions of the asymptotic model are fairly close 
to experiments. When the size of the voids becomes large compared to the distance between 
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Figure 6: Comparison between experimental (A) and theoretical (B) results for quasi-static crack trajectories 
in the case of high value of the ratio between the void diameter and the distance from the void centre to the crack 
trajectory. The numerical simulation correctly predicts the crack path, while the asymptotic approach correctly 
captures only the initial crack deflection. 


the unperturbed crack trajectory and the center of the voids, the approximation error of the 
asymptotic model also increases. Moreover, when the crack propagation reaches a critical 
dynamic regime, kinking of the crack path occurs, and consequently, the crack may be arrested 
as a result of hitting one of the voids. Figs. and illustrate these special cases when the 
crack deviates from the asymptotic model. The experiments have demonstrated that, when 
the underlying hypotheses remain valid, the asymptotic model gives excellent results for a wide 
range of loading rate conditions, as shown in Figs. ii and The importance of this analysis 
lies in the possibility of accurately predicting the failure patterns of brittle materials, such as 
ceramics, with inhomogeneities in the form of small inclusions or voids. 
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